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ABSTRACT 

A multivariate analysis of gamma-ray burst (GRB) bulk properties is presented to 
discriminate between distinct classes of GRBs. Several variables representing burst duration, 
fluence and spectral hardness are considered. Two multivariate clustering procedures are used 
on a sample of 797 bursts from the Third BATSE Catalog: a nonparametric average linkage 
hierarchical agglomerative clustering procedure validated with Wilks' A* and other MANOVA 
tests; and a parametric maximum likelihood model-based clustering procedure assuming 
multinormal populations calculated with the EM Algorithm and validated with the Bayesian 
Information Criterion. 

The two methods yield very similar results. The BATSE GRB population consists of 
three classes with the following Duration/Fluence/Spectrum bulk properties: Class I with 
long/bright/intermediate bursts, Class II with short /hard/ faint bursts, and Class III with 
intermediate/intermediate/soft bursts. One outlier with poor data is also present. Classes I and 
II correspond to those reported by Kouveliotou et al. (1993), but Class III is clearly defined 
here for the first time. 



Subject headings: gamma rays: bursts; methods: statistical; methods: data analysis 
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1. Introduction 

As very few gamma-ray burst (GRB) sources have astronomical counterparts at other wavebands, 
empirical studies of GRBs have been largely restricted to the analysis of their gamma ray properties: 
bulk properties such as fluence and spectral hardness, and evolution of these properties within a burst 
event (Fishman & Meegan 1995). While bursts exhibit a vast range of complex temporal behaviors, their 
bulk properties appear simpler and amenable to straightforward statistical analyses. Studies fall into 
two categories: examination of whether GRB bulk properties comprise a homogeneous population or are 
divided into distinct classes; and search for relationships between bulk properties. Both types of study 
may lead to astrophysical insight, just as the distinction between main sequence stars and red giants and 
the measurement of a luminosity-mass relation along the main sequence assisted the development of stellar 
astrophysics early in the century. 

The most widely accepted taxonomy of GRBs is the division between short-hard and long-soft bursts 
proposed by Dezalay et al. (1992) and Kouveliotou et al. (1993, henceforth K93). K93 noticed a bimodality 
in the burst duration variable Tgg (time within which 90% of the flux arrived), suggesting the presence of two 
distinct types of bursts separated at Tgo — 2 sec. The short bursts have systematically harder gamma-ray 
spectra than longer bursts. The two groups seemed indistinguishable in most other bulk properties, 
although the larger group of long-soft bursts may have a subclass with a different fluence distribution (i.e., 
different < V/Vm >; Katz & Canel 1996) and the groups may have different Galactic latitude distributions 
(Belli 1997). Other researchers point to small groups of bursts with distinctive properties such as the 
soft-gamma repeaters (Norris et al. 1991), two possible classes with differing short-timescale variability 
(Lamb, Graziani & Smith 1993), fast-rise exponential-decay bursts (Bhat et al. 1994), and two types of 
bursts with different ratios of total fluence and >300 keV fluence (Pendleton et al. 1997). 

A variety of relationships between burst properties have also been reported. Norris et al. (1995) find 
an anti-correlation between Tgo (calculated after wavelet thresholding) and peak intensity, consistent with a 
cosmological time dilation. However, a positive correlation between Tgo and total fluence is also seen which 
does not agree with the simplest cosmological interpretation (Lee & Petrosian 1997). Additional reported 

relationships include: Tgo correlated with peak heights (Lestrade 1994), peak energy correlated with peak 
flux (Mallozzi et al. 1995), and peak duration anticorrelatcd with gamma-ray energy (Fenimore et al. 1995). 

Most of these studies suffer from a failure to treat all of the bulk property variables in an unbiased 
and quantitative way. Astronomers typically examine univariate or bivariate distributions, sometimes 
constructing composite variables (such as hardness ratios) with pre-determined relationships to include one 
or two additional variables. But it is quite possible that the complex astrophysics producing GRBs will not 
manifest themselves in simple bivariate plots, just as the division between short-hard and long-soft bursts 
is not evident in spectral variables alone (Pendleton et al. 1994). GRB catalogs, like most multiwavelength 
astronomical catalogs, are multivariate databases and should be treated with multivariate statistical 
methods that can objectively and effectively uncover structure involving many variables (Feigelson & Babu 
1997). Two previous studies take a fully multivariate approach to understanding GRB bulk properties. 
Baumgart (1994) constructs a neural network taxonomy of 99 GRBs from the PVO satellite using 26 
variables representing both bulk burst properties and detailed temporal characteristics (e.g. number of 
peaks, fractal dimension, wavelet transform crossings) and finds two or three distinct GRB classes. Bagoly 
et al. (1997) perform principal components and factor analyses of nine bulk property variables using 625 
GRBs from the BATSE 3B catalog. They find that the relationships in the database are determined 
principally by only three variables: an appropriately weighted fluence, a weighted burst duration, and (to a 



- 3 - 



lesser extent) flux in the highest energy bin. 

We note, however, that it can be dangerous to look for correlations prior to classifying (or establishing 
the homogeneity of) the population. While the anticorrclation between hardness ratio and burst 
duration seen in full samples (K93) may be the manifestation of a single astrophysical process, it may 
alternatively reflect diff'erences between distinct processes. The latter possibility is suggested by a reported 
hardness-duration positive correlation within the long-soft class of bursts (Dezalay et al. 1996; Horack & 
Hakkila 1997). Most multivariate analyses thus begin with a study of homogeneity and classification, and 
then investigate the variance-covariance structure (i.e. correlations) within each class. 

This paper describes a multivariate analysis of GRBs from the Third BATSE Catalog (Meegan et al. 
1996). After defining the sample (§2), we start with a simple statistical description of the variables and 
their bivariate relationships for the entire dataset (§3). We then seek distinct types of clusters in two ways. 
First, a standard nonparametric agglomerative hierarchical clustering analysis is performed (§4) which 
reveals three distinct classes. The statistical significance of the third cluster is validated, under Gaussian 
assumptions, with MANOVA tests. Second, a parametric maximum likelihood model-based clustering 
procedure is adopted which reveals the same three groups and indicates strong evidence for the presence for 
three rather than two groups (§5). The variance-covariance structure of each group is then examined (§6). 
Results are synthesized in the discussion (§7). 

Throughout the paper, we discuss our mathematical techniques to help the reader understand the 
complexities of multivariate analysis. From the vast literature in this subject, we recommend the following 
monographs for interested readers: .Johnson & Wichern (1992) and Jobson (1992) for overviews of applied 
multivariate analysis; Hartigan (1975), Jain & Dubes (1988) and Kaufman & Rousseeuw (1990) for 
multivariate clustering algorithms; Murtagh & Heck (1987) and, more briefly, Babu & Feigelson (1996) and 
Feigelson & Babu (1997), for multivariate methodology in astronomy. 

2. The GRB sample and statistical software 

Our sample is drawn from the Third Catalog of the Burst and Transient Source Experiment (BATSE) 
on board the Compton Gamma Ray Observatory. This 3B catalog has 1122 GRBs detected by BATSE 
between 1991 April 19 and 1994 September 19. The catalog is presented and fully described by Meegan et al. 
(1996). Our database was extracted from the on-line database www.batse.msfc.nasa.gov/data/grb/catalog 
in May 1996, which provides many properties of each burst. There are roughly eleven variables of potential 
astrophysical interest: two measures of location in Galactic coordinates, I and b; two measures of burst 
durations, the times within which 50% (T50) and 90% (Tgo) of the flux arrives; three peak fluxes Pe4, P256 
and P1024 measured in 64 ms, 256 ms and 1024 ms bins respectively; and four time-integrated fluences 
F1-F4 in the 20-50 keV, 50-100 keV, 100-300 keV and > 300 keV spectral channels respectively. Researchers 
commonly consider three composite variables: the total fluence, Ft = Fi + F2 + F3 + F4, and two measures 
of spectral hardness derived from the ratios of channel fluences, H32 = F3/F2 and ^321 = Fs/{Fi + F2). 
Due to the limitations of available multivariate statistical techniques, we ignore other variables of potential 
relevance including the heteroscedastic measurement errors of each quantity (i.e., errors that differ from 
point to point) and truncation values associated with BATSE triggering operations. 

Of the 1122 listed bursts, 807 have data on all the variables described above. Ten bursts listed with 
zero fluences were eliminated. Our sample thus has 797 BATSE GRBs. For some analyses, we also used 
a subset of 644 bursts with 'debiased' durations, T^q. Here the durations are modified to account for the 
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efFect that brighter bursts wiU have signal above the noise for fonger periods than fainter bursts with the 
same time profiles (J. Norris, private communication). 

Statistical analyses in §§3, 4 and 6 were conducted within the Statistical Analysis System SAS/STA'M, 
a very large and widely distributed commercial statistical software package (SAS Institute Inc. 1989). 
S AS/ STAT wocedmes CLUSTER, GLM, PRINCOMP and VARCLUS were used. The analysis in §5 was 
performed with the MCLf/^T software (Banfield & Raftery 1993; Fraley 1998), which is interfaced to the 
Splus statistical package (Splus Version 3.4; MathSoft, Inc. 1996) and its extensions. Further information is 
provided at http://stat.washington.edu/fraley/software.htmj and http: / /lib. stat.cmu.edu/general / mclust . 
For multivariate data visualization, we used the XGobi (Swayne, Cook & Buja 1991) program, available 
from the StatLib software archive at http:/ /lib. stat.cmu.edu/general/XGobi. Hypertext links to a variety of 
public domain software for multivariate analysis, classification and visualization are available at the Penn 



State Stoi Codes Web site, http://www.astro.psu.edu/statcodes 



3. Statistical properties of the entire sample 



We are faced with a multivariate database of 797 objects and 15 variables (11 variables from the catalog, 
3 composite variables, and Tqq). Two initial problems are frequently faced in analyses of multivariate 
databases. First, variables with incompatible units and ranges must be compared. Units can be removed 
by normalization (e.g. replacing Fi by Fi/Ftot), by standardization (e.g. replacing Fi by Fi/a where a is 
the sample standard deviation), or by taking logarithms. Second, the dimensionality of the problem should 
be reduced, as many of the variables are closely interrelated either by construction or by astronomical 
circumstance. Although there are no mathematical rules regulating reduction of dimensionality, it can 
usefully be guided by a correlation matrix showing bivariate relationships and a principal components 
analysis showing multivariate relationships that are mainly responsible for structured variance in the data. 
Scientific reasoning can also be used to eliminate consideration of variables. We conducted a preliminary 
examination of data representations, correlation matrices and bivariate plots, and principal components 
analyses to facilitate choice of variables. When no mathematical preference arose, we selected variables 
most commonly used by previous researchers to facilitate comparison of results. 

Our choices were as follows. We use log variables, rather than normalized or standardized variables. 
We kept information on burst fluence and spectra through Ftot and hardness ratios rather than through 
the original fluences F1-F4. We initially eliminated Pq^ and P1024 from consideration, and later eliminated 
-P256 when we found that its main contribution to the clustering process was noise. We chose to remove 
the location variables {l,b), already established by other researchers to be random for the entire sample, 
but use them later to test for isotropy of subsamples. The debiased Tgo^d is used only in special tests. Our 
analysis was thus performed in six or fewer dimensions using log T50, log Tgo, log Ftot, log ^2567 log H321 
and log H32- 

Tables 1 and 2 give basic statistics for these six variables: means, standard deviations, and bivariate 
values of Pearson's linear correlation coefficient r. For N — 797 and assuming bivariate normal populations, 
any Irl > 0.013 implies that a correlation between the two variables exists at a two-tailed significance level 
P < 0.001 (Beyer 1968, pp. 389 and 283). But from an astrophysical perspective, we might consider any 
relationship with |r| < 0.1 to be of little interest. Figure 1 shows the bivariate scatter plots. 



^ SAS/STAT is a registered trademark of the SAS Institute Inc. 
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The correlation structure of the entire sample (Table 2) shows that the two measures of duration 
and the two measures of spectral hardness have correlation near unity, indicating that they are nearly 
redundant. Ft^i and P256 are quite dissimilar: Fjoj shows a strong correlation with burst duration (e.g. 
Lee & Petrosian 1997) and no relation to hardness, while P256 shows no relation to duration but is mildly 
correlated with hardness (Mallozzi et al. 1995). Burst duration is anticorrelated with hardness (K93; 
Fenimore et al. 1995). The cosmological anticorrelation between duration and peak flux reported by Norris 
et al. (1995) is statistically significant, but accounts for only a few percent of the variance between these 
variables. The correlation matrix based on the debiased Tgo values yields very similar results. 

However, the scatter plots (Figure 1) show a more complex story. First, many plots show inhomogeneous 
distributions inconsistent with the unimodal multinormal (i.e. multivariate Gaussian) population assumed 
by Pearson's r. The distributions often seem bimodal with asymmetrical non-Gaussian shapes. One outlier 
burst is also seen in several projections. We therefore consider the hypothesis that the sample consists of 
two or more distinct classes, and proceed to find the 'clusters' using well-established methods. 

4. Nonparametric hierarchical cluster analysis 

4.1. Methodological background 

Agglomcrative hierarchical clustering is a procedure based on the successive merging of proximate pairs 
of clusters of objects. It produces a clustering tree or dendrogram starting with N clusters of 1 member (or 
a coarse partition based on prior knowledge) and ending with one cluster of N members. Unfortunately, 
there are many possible ways to proceed; mathematics provides little guidance among the choices and no 
probabilistic evaluation of the results without the imposition of additional assumptions. The scientist must 
make four decisions to fully define the clustering procedure: 

1. Creating unit-free variables is essential for meaningful treatment of objects in multivariate space 
(§3). A favorite choice by statisticians is standardization, where each variable is normalized by the 
standard deviation of the sample. Astronomers more commonly make logarithmic transformations or 
construct ratios of variables sharing the same units. We follow the tradition of ORB researchers by 
measuring spectral hardness with ratios of fluences having the same units, and making logarithmic 
transformations of all variables. 

2. The metric defines the meaning of proximity between two objects or clusters. Common choices are 
the simple Euclidean distance between unit-free variables and the squares of Euclidean distances. We 
chose the former option for most of the analyses in this section. 

3. Several merging procedures can be used. One might begin by merging the clusters with the nearest 
neighbors. This is called Single Linkage clustering and is most familiar in astronomy where it is 
frequently called the friends-of-friends algorithm. It tends to produce long stringy clusters, and is 
equivalent to a well-known divisive clustering procedure known as pruning the minimal spanning 
tree. Complete Linkage proceeds by maximizing the distance between clusters, and leads to evenly 
bifurcating dendrograms. For most of our analysis, we choose Average Linkage where the distance 
between two clusters is the average of the distances between pairs of observations, where each member 
of the pair comes from a different cluster. This is a compromise between Single and Complete Linkage 
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and tends to give compact clusters. Specifically, the distance between clusters K and L is given by 
(e.g. SAS Institute Inc. 1989, pp. 529ff; Johnson & Wichern 1992, pp. 584ff.) 



-DifL = |xk - XL H \ , (1) 

UK riL 

where the bar indicates an unweighted mean, Wk = l^j — x|^, and Uk is the number of members 

of the fc-th cluster. Another popular choice is Ward's minimum variance criterion where the distance 
between the two clusters is the ANOVA (analysis of variance) sum of squares between two clusters 
added up over all variables (Ward 1963), 

Dkl = \^k-^l\V{— + —)- (2) 

UK riL 

If the sample is generated by a mixture of multinormal (i.e. multidimensional Gaussian) distributions 
where each distribution has covariance matrix of the form S^/, this method joins clusters to maximize 
the likelihood at each level of the hierarchy, and so is a special case of the model-based clustering 
methodology to be discussed in §5. 

As the procedure gives a hierarchy from N clusters with 1 object down to one cluster with N objects, 

the user must choose how many clusters to report as scientifically important clusters. This choice 
can be assisted by examination of two statistics. The squared correlation coefficient, i?^, states the 
fraction of the total variance accounted for by a partition into g clusters. 



^^ = l- ^iv „ (3) 



j=il^i - ^1 

The squared semi-partial correlation coefficient, measures the difference in the variance between 
the resulting cluster and the immediate parent clusters normalized by the total sample variance, 

„2 _ Wm-Wk-Wl 
E,^i|xi-x|2 • 

i?^ thus tells how much of the scatter is explained by a given level of clustering, and R^p tells how 
much improvement is achieved between levels. 

We emphasize again that there is no mathematically 'best' choice, although extensive experience with 
problems in many fields has led to a preference for certain combinations (e.g. standardized variables and 
Ward's minimum variance criterion). We conducted extensive experiments with different choices. 



4.2. Results 



The last several levels of the clustering tree for the 797 GRBs using the six unit-free variables shown in 
Table 1, Average Linkage and a Euclidean metric are shown in Figure 2 (left panel) with details in Table 3a. 
The action taken at each level is indicated in column 2 of Table 3, which may refer to a level higher in the 
tree which (for brevity) is not shown here. Two types of mergers arc seen: the incorporation of 'twigs' of 
one or a few GRBs into a large preexisting 'trunk' (levels 1, 3, 4, 5 and 7); and the union of two substantial 
branches into a single larger trunk (levels 2, 6 and 8). The first type has little effect on the variance of the 
sample with R^ < 1%. The single GRB brought into the main trunk at level 1 is the distant outlier seen 
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in several panels of Figure 1. The level 2 merger of clusters with 190 and 606 members is clearly the most 
important structure, accounting for roughly 53% of the variance of the entire sample. This is the bifurcation 
of the sample into two classes easily seen in Figure 1 and noted by K93 and others. The principal finding 
that is not immediately obvious from Figure 1 is the structure indicated at level 6. The main trunk of 599 

bursts (plus a few twigs to be merged later) is divided into groups of 93 and 506 bursts. This division 
accounts for 10% of the total variance of the sample, indicated in both the and R^p values. 

We found that the twigs in the tree structure disappear if the peak flux -P256 variable is omitted and 
the analysis is made in 5-dimensional space (Figure 2, right panel, and Table 3b). Here the largest cluster 
of 593 members is formed by the union of clusters with 107 and 486 bursts, again accounting for 10% of the 
sample variance. It is possible that P256 is a nuisance variable irrelevant to the basic astrophysics of GRBs, 
producing noisy 'twigs' seen in Table 3a and Figure 2 (left panel). 

We tested many variants of hierarchical clustering. We replaced Average Linkage hypothesis 
with Complete Linkage, Single Linkage and Ward's minimum variance criterion. The Ward's criterion 
computation, for example, gave three clusters with 468, 145 and 184 bursts. We clustered using 
nonparametric density estimation based on the 100 nearest neighbors, and clustered using the principal 
components rather than the observed variables. Various methods were tried with both the observed T90 
values and debiased Tgg values, with little effect on the results. All methods showed two strong clusters and 
the outlier but, in some cases, the third cluster appeared only weakly. 

To proceed further, we choose a single clustering structure for detailed study: the 5-dimensional 
Average Linkage analysis (Table 3b) with three clusters — Class I with 486 bursts, Class II with 203 
bursts. Class III with 107 bursts. Class IV consisting of the single outlier is ignored because of independent 
evidence that its properties are due to data of poor quality (§7). The membership of these clusters is given 
in Table 4, and four projections of the clusters onto two-dimensional scatter plots are shown in Figure 3. 
These are frames from the 'grand tour' movie of the 5-dimensional dataset provided by the XGobi software 
where each cluster is 'brushed' with a different symbol. Note that, in general, there is no reason why 
classification structure should be most evident in projections parallel to the variable axes shown in Figure 
1. It is more important that the clusters show cohesion in many projections of the dataset. The grand 
tour of the 797 GRBs shows that Classes I, II and the outlier are very distinct in most projections. Class 
III often lies between Classes I and II (e.g. top panels of Figure 3), but in other projections is offset from 
the line between Classes I and II (e.g. bottom panels of Figure 3). It also appears elongated along some 




Fig. 2. — Diagram of the base of the dendrogram of Average Linkage hierarchical clustering procedures in 6 
dimensions (left panel) and 5 dimensions (right panel). The number of members in each branch is indicated 
(see Table 3). Class IV is the spurious outlier. 
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Figure 3. Four snapshots from the XGobi grand tour of the 5- dimensional database 
with bursts brushed according to the nonparametric Average Linkage clustering re- 
sults: Class I (•), Class II (x), Class III (□) and Class IV (O)- Panel (a, upper right) 
shows the projections of the five axes, which are suppressed in the other panels. 
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projections, while the larger Classes I and II appear roughly hyperspherical. 

This analyses described here provide considerable evidence for three major clusters and an outlier. But 
as some nonparametric clustering procedures did not find a strong third cluster, there is some some worry 
that Class III is simply a group of bursts with properties intermediate between the Classes I and II. While 
nonparametric hierarchical clustering methods can not address this question, it can be investigated with 
parametric methods. 

4.3. Validation of the classification 

Mathematically well-founded methods for evaluating the statistical significance of a proposed 
multivariate classification scheme are available under the assumption that the population is a multinormal 
mixture; that is, the objects of each class are drawn from multivariate Gaussians. All relationships 
between the variables must thus be linear (as in Table 2), although the relationships may differ between 
clusters. There is no requirement of sphericity, so that clusters may have shapes akin to pancakes or cigars 
with arbitrary orientations in multidimensional space. The separate existence of each of the postulated 
subpopulations can be tested using multivariate analysis of variance (MANOVA). 

The model can be expressed as follows (e.g. Johnson & Winchern 1992, pp. 246ff). For a p-dimensional 
dataset of g clusters each with n; members, the i-th GRB in the j-th. cluster gives a p-dimensional vector 

Xjj = n + Tj + Cij (5) 

where ii is the overall population mean, Tj is the offset of the j-ih cluster mean from /k, and Cij are 
independent normal variables with zero mean representing the scatter of individual points about the mean. 
We test the null hypothesis 

Ho:ti=T2 = ... = Tg = (6) 

that the cluster means are not offset from each other. We construct two matrices of sums of squares and 
cross-products as follows: 

9 

B = ^n;(x; -x)(x; -x)"^ 

1=1 

g ni 

W = 5^5^(xy-Xi)(xy-Xif, (7) 

1=1 3 = 1 

where ^ is the vector transpose. Three test statistics have been proposed to test the null hypothesis (e.g. 
SAS Institute Inc. 1989, pp. 17ff)): 

Wilks' Lambda A* = det(W)/dct(B + W), 

Pillai's trace V = trace[B(B + W)"^], and (8) 
Hotelling — Lawley's trace U = trace(W~^B). 

The distributions of these statistics have been determined mathematically. For example, for large 

= X]f=i t;he quantity —{N — 1 — [p + 5)/2)lnA* has approximately a chi-squared distribution with 
p{g — 1) degrees of freedom (Wilks 1932; Bartlett 1938). More generally, the distributions are related to 
the non-central F distribution. For the 2-sample case, Hotelling-Lawley's trace is commonly known as the 
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Mahalanobis statistic. One can thus accept or reject the null hypothesis that the clusters have the same 
mean location at a chosen level of statistical significance. 

The results of our MANOVA calculations are summarized in Table 5. The columns give the values 
of the three MANOVA statistics, followed by details for the Wilks' A*: the corresponding value of the F 
statistic, the numerator and denominator degrees of freedom for that F value; and the resulting P value. 
Details for Pillai's and Hotelling-Lawley's traces are omitted, but give similar results in all cases. The first 
row tests the null hypothesis that the Classes I, II and III have the same mean, the second row tests the 
equality of Classes I and II, and so forth. The F values are very high in all rases, indicating that the 
clusters are different with extremely high statistical significance {P « 10^^). □ This is a clear quantitative 
demonstration that at least two clusters exist among GRBs (for a univariate test, see Ashman, Bird & Zepf 
1994), which was qualitatively reported by Delazay et al. (1992) and K93. The other rows in Table 5 test 
the hypotheses that each proposed class has the same mean as each other class. 

One problem with these MANOVA tests is that they are conditional on the classification that has been 
found using the clustering algorithm. Because the clustering algorithm is constructed to find groups that 
are different from one another, tests such as these tend to be biased towards finding structure, perhaps 
where none exists. Although the MANOVA results seem to indicate very strong evidence of structure, they 
can not be taken as definitive for this reason. Tests arising from model-based clustering can overcome this 
problem, as discussed in the following section. 

5. Model-based m£Lximum likelihood clustering analysis 

5.1. Methodological background 

In the previous section, we conducted a hierarchical clustering analysis without making assumptions 
regarding the shapes of the clusters, but needed the parametric assumption of normality to estimate 
the statistical significance of the resulting classification. It is reasonable to conduct the entire analysis, 
both clustering and validation, within a model-based framework. We report here an analysis of this 
type again assuming that the GRB population consists of a mixture of multivariate Gaussian classes. 
Early development of this model for clustering is discussed in McLachlan & Basford (1988); we use more 
recent developments here. First, an initial classification for each possible number of clusters is found via 
agglomerative hierarchical clustering (Murtagh & Raftery 1984, Banfield & Raftery 1993, Fraley 1998). 
Next, the EM (Expectation-Maximization) algorithm is used to refine partitions obtained from hierarchical 
clustering (Celeux & Govaert 1995, Dasgupta & Raftery 1998). Finally, the Bayesian Information Criterion 
is used to select the 'best' partition among those associated with different numbers of clusters (Dasgupta & 
Raftery 1998). 

In the model considered here, the p-dimensional observations are drawn from g multinormal groups, 
each of which is characterized by a vector of parameters Ok for k = 1, . . . ,g. Our goals are: to determine 
the number of GRB types, g; to determine the cluster assignment of each burst; and to estimate the mean 
Ilk and covariance matrix for each cluster. Following Fraley (1998), the density of an observation 



^Note that it is not meaningful to quote probabilities like P = 10 * as the tails of the distribution are poorly determined 
unless the sample size is extremely large. 
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from the A;th subpopulation is expressed as follows: 

M^i\Ok) ^ MVN{^,k,'Sk) k = l,...,g (9) 

where MVN means multivariate normal. We estimate the parameters using the principle of maximum 
likelihood. In the hierarchical clustering phase, we use the classification likelihood 

Lo{e,i\^) = Jlf^,{^i\6^,) (10) 

i=l 

where x = Xi,X2, . . . ,XAr represents the observations and 7 = {71,72, . . . ,7jv} is the cluster assignment: 
"fi — k when Xi comes from the k-th group. Equivalcntly, 

Lc{ni, . . . , /ife; Si, . . . , Sfe I x) = 

n n (27rr/2|Sfe|-V2 exp{-l(x, - /XferS^i(x, - /.fc)}, (11) 

where /fe = {i : 7, = k} is the set of indices corresponding to observations belonging to the fc-th group. 

The method used here for maximizing the likelihood (Fraley 1998) and implemented in the MCLUST 
code involves parameterization of the Sfe matrices in terms of their eigenvectors and eigenvalues (analogous 
to a principal components analysis), and iterative relocation of the clusters using the EM Algorithm. The 
EM (Estimation-Maximization) Algorithm (Dempster, Laird & Rubin 1977), one of the most successful 
methods in modern statistics, is a procedure for itcrativcly maximizing likelihoods in a wide variety of 
circumstances. For example, the Lucy-Richardson algorithm in astronomical image restoration is the EM 
Algorithm. In the present application, we apply EM to the mixture likelihood 

N g g 

L„(^,7|x) = n$^T-fe/fe(xi|^fc); 5^rfe = l, (12) 

i=l fc=l fe=l 

where are mixing probabilities associated with each group. For a given number g of components in 
the mixture, we use EM to estimate the conditional probability that observation x^ belongs to the k-th 
group for each i and selected k via maximum likelihood. Although the computational procedure has some 
limitations (e.g. convergence of the EM iterations is not guaranteed; clusters cannot be extremely small), it 
is generally efficient and effective for Gaussian clustering problems when started from reasonable partitions 
such as those produced by hierarchical agglomeration. 

We use the Bayes factor to assess the evidence for a given number of clusters against a different number 
of clusters. The Bayes factor, defined in the context of Bayesian statistics, is the posterior odds for one 
model against the other when the prior odds are equal to one (i.e., when one does not favor one model 
over the other a priori). Kass & Raftery (1995) review the use of Bayes factors in adjudicating between 
competing scientific hypotheses on the basis of data. The Bayes factor for a model M2 against a competing 
model Ml (say, for three vs. two classes of GRBs) is defined as 

p(x|M2) 

Bayes factor = , ' , , , 13 
p(x|Mi)' ^ ' 

where p(x|Mj) for j = 1,2 is obtained by integrating the likelihood times the prior density over the 
parameters of the model. It can be viewed as a likelihood ratio, but it differs from the usual frequentist 
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ratio that underlies the hkehhood ratio test in that the latter is obtained by maximizing (rather than 
integrating) the likelihood over the model parameters. 

Twice the logarithm of the Bayes factor can be approximated by the Bayesian Information Criterion 
or BIG (Schwarz 1978), 

BIG = 2(Zi - h) - (mi - m2)logN, (14) 

where h is the likelihood and mi is the number of parameters for one mixture model, and similarly for I2 

and 7712 ■ The BIG measures the balance between the improvement in the likelihood and the number of 
model parameters needed to achieve that likelihood. While the absolute value of the BIG is not informative, 
differences between the BIG values for two competing models provide estimates of the evidence in the data 
for one model against another. Gonventionally, BIG differences < 2 represent weak evidence, differences 
between 2 and 6 represent positive evidence, 6—10 strong evidence, and > 10 very strong evidence (Jeffreys 
1961, Appendix B; Kass & Raftery 1995). The use of the BIG in choosing clusters in a mixture or clustering 
model is discussed by Roedcr & Wasscrman (1997) and Dasgupta & Raftery (1998). 

Bayes factors and BIG have the advantage that they can be used to assess the evidence for a null 
hypothesis, unlike standard significance tests which can only reject a null hypothesis. They can also easily 
be used to compare non-nested models, again unlike standard significance tests which require competing 
models to be nested. 

5.2. Results and validation 

To reduce the dimensionality of the problem and the complexity of the calculation, we eliminated the 
highly redundant T50 and H32 variables (see Figure 1) and considered only the 3 variables T90, Ftot and 
H321 for the sample of 797 BATSE GRBs. The MCLUST model-based clustering procedure described 
above was run for trials of g = 1, 2, .... 24 groups. The resulting values of BIG(.g) are plotted in Figure 4. 
The maximum BIG is achieved for three classes. Most importantly, the BIG value for (7 = 3 is ~ 68 units 
above that for 3 = 2. This corresponds to strong evidence indeed for the presence of three groups rather 
than two. This result strongly confirms the analysis in §4 indicating the existence of three clusters, and this 
time is free of the problem that the MANOVA tests are conditional on the estimated partition. The result 
here takes account of the fact that the partition is not known in advance. 

We have also calculated the BIG for g = 1, . . . , 9 with various constraints on the covariance matrix S 
such as hypersphericity and uniformly shaped ellipsoids. Spherical clusters give poor fits. Uniform ellipsoids 
give good fits with 4 and 8 clusters. But in all cases, the maximum likelihood assuming 2 clusters is much 
lower than the likelihood of > 3 clusters. 

The cluster assignment vector 7 for the 5 = 3 model with unconstrained E is given in Table 4. Over 
85% of the assignments are the same as those obtained from the nonparametric hierarchical clustering 
procedure in §4, so that we note only differences between the two clustering results with * and t markings. 
All but one of the 96 assignment differences move bursts from Glasses I and II into Glass III. The close 
agreement between the cluster assignments in the two methods reinforces confidences in the conclusions 
from both of them. 
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Fig. 4. — Bayesian Information Criterion value vs. assumed number of clusters for the model-based 
Gaussian mixture maximum likelihood clustering analysis 
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Fig. 5. — Burst types and Tgo burst duration based on 5-dimensional Average Linkage clustering 
procedure. Class I = white, Class II = grey and Class III = black. The Class IV outlier is shown as 
an X. 
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6. Cluster properties 

We can now examine the properties of GRBs within each cluster with reasonable confidence that the 
populations are distinct from each other but internally homogeneous. These properties become inputs to 
astrophysical theories seeking to explain GRB bulk properties. Table 6a lists the means and standard 
deviations of the principal variables for each cluster based on both the nonparametric and model-based 
clustering procedures. The two methods give very similar results. The three types are well-separated in the 
burst duration variables: Cluster I bursts have the longest durations around 10-20 seconds, Cluster II bursts 
have the shortest durations below 1 second, and Cluster III bursts have intermediate durations around 2 
seconds. This is shown clearly in Figure 5, which projects each class onto the univariate Tgo axis. Cluster 
III bursts are also intermediate in their fluences, although their fluence distribution overlaps that of the 
fainter Class II bursts. The hardness ratios of all three clusters overlap considerably, but Class III bursts 
have the softest spectra and Class I bursts have intermediate spectra. We can thus classify the types in 
the three principal dimensions Duration/Fluence/Spectrum (Table 6c): Class I is long/bright/intermediate. 
Class II is short/faint/hard, and Class III is intermediate/intermediate/soft. 

A major constraint for the astrophysical interpretation of GRBs has been the remarkable isotropy of 
their spatial distribution in the celestial sphere. It is possible that, while the bulk of GRBs are isotropic 
and have an inferred extragalactic origin, some class of GRBs have significant anisotropy which would 
reflect a Galactic origin (see Lamb 1995). We apply four statistical tests for isotropy discussed by Briggs 
(1993) and applied by Briggs et al. (1996) to various subsamples of the 3rd BATSE Catalog of GRBs. The 
statistics are: <cos^ >, where 9 is the angle between a burst and the Galactic center; <sin^6 ~ I where 
b is the Galactic latitude; Raylcigh- Watson W; and Bingham B. <cos9 > tests the dipole moment around 
the Galactic center, <sin^6 — | > tests the quadrapole moment with respect to the Galactic plane, W tests 
the dipole moment around any point in the celestial sphere, and B tests the quadrapole moment around 
any plane or two poles. The expected values for the four statistics assuming random isotropic distribution 
on the sphere are 0, 0, 3 and 5 respectively. The asymptotic distributions of these statistics are known. 

Table 6b shows the results of this analysis for Clusters I-III, kindly calculated for us by Michael Briggs. 
No deviations from isotropy are found. The <cos9 > and <sin^ ~ I > values lie within one standard 
deviation of the expected value for a random distribution. The W and B values must be larger than the 
expected value to indicate anisotropy. The only such case, Class II with B =7.32, has a deviation with 
very low significance {Prob < 0.2). We thus do not confirm Belli's (1997) report of significant differences in 
spatial distributions of burst Classes I and III, although we did not specifically test the Galactic latitude 
distribution. 

In principle, the relative populations of the three classes may be an important constraint on 
astrophysical theory. We find that Class I contains more than half of the bursts with the remainder divided 

between Class II and Class III (Table 6c). But wc do not believe our analysis gives a precise census for two 
reasons. First, the exact assignments of individual bursts to clusters depend on the detailed assumptions of 
the clustering algorithms. For example. Class II is larger than Class III in the 5-dimensional nonparametric 
procedure but is smaller in the 3-dimensional model-based procedure. Second, the numbers of weaker bursts 
in Classes II and III are strongly dependent on the details of the BATSE instrument's burst triggering 
process which produces a complicated tnmcation bias for fainter bursts. 

We look for structure within each of the clusters by computing correlation coefficients similar to those 
in §3 for the entire sample. Results are given in Table 7. Here we see a systematic difference between the 
two clustering methodologies: nonparametric Average Linkage clustering tends to give stronger correlations 
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between the variables than the model-based clustering. For example, in the nonparametric analysis we find 
significant positive correlations between total fluence and hardness in Classes I and II, and a correlation 
between duration and fluence in Class II. However, we attach more credence to the model-based results for 
this purpose than the Average Linkage results because the former method is specifically designed to provide 
optimal estimates of the within-group covariances given the clustering model. The model-based results do 
not give strong evidence for any non-zero correlations between variables, suggesting that the partition into 
three clusters explains all of the correlation between variables in the full dataset. 

7. Discussion 

We thus find, using clustering and validation methods with different mathematical underpinnings, that 
three classes of GRBs are present in our large subset of the Third BATSE Catalog. Most of the structure 
can be found using three fundamental burst properties, Duration/Fluence/Spectrum. The class properties 
and relation to previous research can be briefly summarized as follows: 

Class I These long/bright/intermediate bursts correspond to the well-known populous long-soft class of 
K93 and others. Within this group, we do not confirm a hardness-duration correlation reported by 
Dezalay et al. (1996) and Horack & Hakkila (1997). 

Class II This short/faint/hard group corresponds to the short-hard burst type of K93 and others. 

Fluence-duration and fluence-hardness correlations may tentatively be present within the class. Note 
that while the mean location of this type is consistent in the two clustering schemes, its size and 
population (e.g. 1/2 or 1/4 that of Class I) differs between clustering algorithms. 

Class III The discovery of this group with intermediate/intermediate/soft properties is the principal 

result of this study. The group is easily distinguished in the projections of Figure 3, but can also be 
discerned in some panels of Figure 1. For example, it lies between Classes I and II in the T^q — H32, 
T90 — Ftot and Tgo — scatter plots. In the univariate T90 distribution shown in Figure 5, Class 
III accounts for most, but not all, of the bursts in the small peak around 2 < Tgo < 5 sec between 
the major short and long duration peaks. It is possible that our Class III is related to the class of 
no-high-energy (NHE) burst and peaks discussed by Pendleton et al. (1997). These bursts have 
unusually weak F4 emission, soft 50 — 300 keV spectra, and low Ftot- However, the NHE class does 
not appear to exhibit a clear duration segregation from other bursts as we find for Class HI. Class 
HI does not appear to be the third cluster found by Baumgart (1994, see his Table 3), but the high 
dimensionality of his analysis prevents a simple comparison with our low dimensionality study. 

Outlier BATSE trigger event 2757, burst 3B 940114, is the outlier in the nonparametric analysis of §4 and 
is clearly visible in many projections in Figures 1 and 3.i It has an exceedingly soft hardness ratio 
and short burst duration. But examination of the original BATSE database shows that the Fi — F4 
fluxes are very weak with large measurement uncertainties. The published 3B catalog gives only 
an upper limit to its total fluence and no estimate of its hardness ratio (Meegan et al. 1996). The 
unusual properties of this burst are thus illusory and are due to its very weak fluence. 



^The model-based analysis of §5 cannot locate clusters with very few members and assigned this event to Class II. An 
extension of model-based clustering that models outliers as Poisson noise can do this (see Banfield and Raftery 1993, Dasgupta 
and Raftery 1998), but it does not seem necessary in this application. 
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The multivariate analysis described here is not comprehensive and may not have uncovered all of 
the structure in the Third BATSE Catalog of bulk GRB properties. Our reduction of dimensionality 
may have been too severe, omitting, for example, the potentially important as a distinct variable 
(Pendleton et al. 1997; Bagoly et al. 1997). Many methodological options were not exercised. For example, 
it would be valuable to repeatedly apply the /c-means partitioning algorithm to the database under the 
assumption that three clusters are present (see Murtagh 1992 for an astronomical application of this 
method), check for skewness or kurtosis in the clusters, and undertake an oblique decision tree analysis 
to give analytical formulation to hyperplanes separating the clusters (see White 1997). Codes for these 
and many other multivariate techniques are publicly available through the Web metasite StatCodes at 
www.astro.psu.edu/statcodcs. 

However, the efforts described here are far more capable of finding and quantifying clustering in the 
database than most previous analyses (§1). Previous studies have been based on qualitative rather than 
quantitative procedures for identifying structures, and provide no statistical validation of their claims. 
It is thus not surprising that we uncovered structure missed by previous researchers. In particular, our 
confidence in the presence of a third cluster. Class III, is strong. Two completely independent mathematical 
procedures (§4 and §5) found very similar structure, each validated with high statistical confidence. 

It is possible that the clustering reported here is indeed present in the database, but does not have an 
astrophysical origin. The complex triggering mechanism of the BATSE instrument mechanism and biases 

in bulk property values at low signal-to-noise ratio are two problems that probably affect the multivariate 
structure. We have investigated one manifestation of the latter effect using Tqq and found no effect on our 
results. Instrumental biases generally affect the number of bursts found in some regions of the multivariate 
hyperspace (thereby biasing log N - log S distributions) and may alter the location of clusters, but are 
unlikely to cause the appearance of clustering that is not present in the underlying population. Nonetheless, 
since the BATSE instrument identifies bursts on three separate timescales, it is possible that the third 
cluster here is related to a selection effect associated with the BATSE triggering mechanisms. 

We conclude that the Third BATSE Catalog shows three statistically significant types of bursts 
(Duration/Flucncc/Spectrum): Class I GRBs are long/bright/intermediate. Class II GRBs are 
short/faint/hard, and Class III GRBs are intcrmediate/intcrmediate/soft. Unless the separation of Class III 
from the other types is due to some subtle BATSE instrumental effect, these types are likely to be real and 
their existence should be considered an important input into astrophysical theories for GRBs. For example, 
the three types may reflect different types of external environments and internal shocks in relativistic fireball 
models (Meszaros & Rees 1993; Panaitescu & Meszaros 1998). Note that statistical anlaysis is unable 
to determine whether burst types represent fundamentally different astrophysical processes or distinct 
conditions within a single astrophysical model. 

Our results can be confirmed and extended in two fashions. First, the analysis described here can be 

validated with several hundred more bursts collected by BATSE since the September 1994 cutoff in the 
database used here. Second, following Baumgart (1994), the dimensionality of the problem can be enlarged 
to include detailed characteristics of the burst temporal behaviors. Burst smoothness vs. peakiness, 
characteristic wavelet scales, spectral evolution, and other parameters can be included. With this enlarged 
database, one can perform both an unsupervised exploratory cluster analysis similar to that described 
here, and MANOVA-type analyses that assume the existence of the three groups to determine whether the 
clusters have distinctive temporal properties. 
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Table 1: Average GRB properties for the entire sample 



Variable 


Mean 


S.D" 


log T50 (sec) 


0.55 


0.92 


log T90 (sec) 


0.96 


0.92 


log Ftot (erg cni"2) 


-5.61 


0.76 


log -P256 (photons s~^ cm~^) 


0.16 


0.45 


log ^321 


0.25 


0.33 


log H32 


0.48 


0.30 



" Sample standard deviation 



Table 2: Correlation coefficients for the entire sample 





log ^50 


log T90 


log Ftot 


log P256 


log H32I 


log ^^32 


log T50 


1.00 












log T90 


0.97 


1.00 










log Ftot 


0.63 


0.66 


1.00 








log P256 


-0.01 


0.04 


0.59 


1.00 






log H321 


-0.36 


-0.36 


0.02 


0.24 


1.00 




log H32 


-0.35 


-0.35 


-0.00 


0.19 


0.96 


1.00 



Table 3: Average linkage hierarchical cluster analysis 



Level 


Merger 


Members 


r>2 


i?2 




(a) Six-dimensional an 


alysis 




8 


10 -1- 15 


506 


0.08 


0.65 


7 


14 + 137 


93 


0.00 


0.65 


6 


8 + 7 


599 


0.10 


0.55 


5 


9 + 266 


188 


0.00 


0.55 


4 


5 + 26 


190 


0.00 


0.55 


3 


6 + 12 


606 


0.01 


0.54 


2 


3 4 


796 


0.53 


0.01 


1 


2 + 616 


797 


0.00 


0.00 




(b) Five-dimensional analysis 




6 


15 + 21 


107 


0.01 


0.70 


5 


10 + 8 


486 


0.01 


0.69 


4 


7 + 20 


203 


0.01 


0.68 


3 


6 + 5 


593 


0.10 


0.58 


2 


3 + 4 


796 


0.58 


0.00 


1 


2 + 616 


797 


0.01 


0.00 
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TABLE 4 

Burst classes from Average Linkage clustering (3B Trigger number) 



(a) Class I (486 bursts) 



107* 


472 


816 


1159* 


1533 


1657* 


1982 


2138 


2304 


2431* 


2551 


2703 


2877 


3003 


3109 


109 


473 


820 


1192 


1540 


1660 


1989 


2140 


2306* 


2432 


2560 


2706 


2889 


3005 


3110 


110 


503 


824 


1196 


1541 


1661 


1993 


2143 


2309 


2435 


2569 


2709 


2890 


3011 


3115 


111 


540* 


825 


1197 


1546 


1663 


1997 


2148 


2310 


2436 


2570 


2711 


2891 


3012* 


3119 


114 


543* 


829 


1200 


1551 


1667 


2018 


2149 


2311 


2437* 


2581 


2725* 


2894* 


3015 


3120 


121 


548 


840 


1213 


1552 


1676 


2019 


2151 


2315 


2438 


2586 


2727 


2897 


3017 


3128 


130 


549 


841 


1218* 


1558 


1683* 


2037 


2156 


2316 


2440 


2589 


2728* 


2898 


3026* 


3129 


133 


559 


867* 


1235 


1559 


1687 


2041* 


2181 


2321 


2441 


2593 


2736* 


2900* 


3029 


3130 


143 


563 


869 


1244 


1561 


1700* 


2044* 


2187* 


2324 


2443 


2600 


2749* 


2901* 


3032 


3131* 


148 


577 


907 


1279 


1567 


1704 


2047 


2188 


2325 


2446 


2603 


2751 


2913 


3035 


3132 


160 


591 


927* 


1288 


1574 


1709* 


2053* 


2189 


2328 


2447 


2606 


2753 


2916 


3039* 


3134 


171 


594 


938 


1291 


1578 


1711 


2061 


2190 


2329 


2450 


2608 


2770 


2919 


3040 


3135 


204 


606 


946 


1303 


1579 


1712 


2067* 


2191 


2340 


2451 


2610 


2774 


2922 


3042 


3136* 


211 


630 


973 


1318 


1580 


1714 


2069 


2193 


2344 


2452 


2611 


2775 


2924 


3055 


3138* 


214 


647 


999* 


1384 


1586 


1717 


2070 


2197 


2345 


2472 


2619 


2780 


2925 


3056 


3142 


219 


658 


1009 


1385 


1590 


1730 


2074 


2202* 


2346 


2476 


2620 


2790 


2927 


3057 


3143* 


222 


659 


1025* 


1390 


1601 


1731 


2077 


2203* 


2362 


2477 


2628* 


2793 


2929 


3067 


3153 


223 


660 


1036 


1396 


1604 


1733 


2079 


2204 


2367* 


2482 


2634 


2797 


2931 


3070 


3155* 


226 


673 


1039 


1406 


1606 


1734* 


2080 


2211 


2371 


2484* 


2636 


2798 


2932* 


3071 


3156 


235 


676 


1042 


1419 


1609 


1740 


2081 


2213 


2373 


2495 


2640* 


2799* 


2947 


3072* 


3159 


237 


678 


1046 


1425* 


1611 


1742 


2083 


2219* 


2375 


2496 


2660 


2812 


2948 


3074 


3164 


249 


685* 


1085 


1432 


1614 


1806 


2087* 


2228 


2380 


2500 


2662 


2815 


2950 


3075 


3168 


257 


686* 


1086 


1440 


1623 


1807 


2090 


2232 


2383 


2505 


2663 


2825 


2953 


3076 


3171 


288 


692 


1087 


1446 


1625 


1815 


2093 


2233 


2385* 


2508 


2664 


2831 


2958 


3080 


3174 


332 


704 


1122 


1447 


1626 


1819 


2101 


2230 


2387 


2510 


2665 


2843 


2961 


3084 




351 


717* 


1123 


1449 


1628 


1830 


2102* 


2244 


2391 


2511 


2671* 


2852 


2984 


3085 




394 


741* 


1126 


1452* 


1642 


1883 


2106 


2252 


2392 


2519 


2681 


2853* 


2985 


3091 




398 


761 


1141 


1456 


1646 


1885 


2110 


2253 


2394 


2522 


2688 


2855 


2992 


3093 




404 


764 


1148 


1458 


1651 


1886 


2111 


2267 


2405 


2528 


2691 


2856 


2993 


3100 




408 


773 


1150 


1467 


1652 


1922* 


2112 


2276 


2419 


2530 


2695 


2857 


2994 


3101 




451 


795 


1152 


1468 


1653 


1924 


2119* 


2277 


2428 


2533 


2696 


2862 


2996 


3102 




467 


803* 


1156 


1472 


1655 


1956 


2122 


2287 


2429 


2537 


2697 


2863 


2998 


3103 




469 


815 


1157 


1515 


1656 


1967 


2133 


2298 


2430 


2541 


2700 


2864 


3001 


3105 
















(b) Class 11 (203 bursts) 














138 


512 


856 


1154 


1635 


2003 


2146 


2291 


2384 


2523 


2693 


2846 


2975 


3094 




185 


537* 


878 


1211 


1636* 


2040* 


2155 


2312 


2395 


2529 


2701 


2849 


2977 


3113* 




207 


547 


906 


1223 


1659 


2043 


2159 


2317* 


2434 


2536 


2715* 


2851* 


2978 


3114* 




218* 


551 


909 


2389 


1662 


2049 


2161 


2320 


2448 


2564 


2748 


2860* 


2987 


3118 




229 


568 


936* 


1308 


1665 


2068 


2163 


2326 


2449 


2583 


2755 


2873 


2988* 


3121 




254 


575 


1051 


1359 


1680 


2095 


2167 


2327* 


2454* 


2585 


2788 


2879 


2995 


3137 




289 


603* 


1073 


1404* 


1694 


2099* 


2201 


2330 


2463 


2597 


2795 


2892 


3027 


3152 




297 


677 


1076 


1453 


1719* 


2103 


2205 


2332 


2464 


2599 


2800 


2896 


3037 


3173 




298 


729 


1088 


1461 


1736 


2115 


2206 


2352* 


2485 


2614 


2801 


2910 


3038 






432 


788 


1096 


1463 


1741 


2117 


2217* 


2353 


2487 


2615 


2810* 


2918 


3043 






444 


799 


1097 


1481 


1760 


2125 


2220 


2357 


2502 


2623* 


2814 


2933 


3051 






474 


809 


1102 


1518 


1791 


2126 


2265 


2360 


2504 


2632* 


2821 


2952 


3066 






480 


830 


1112 


1553* 


1851 


2132 


2268* 


2365* 


2512 


2649 


2823 


2964 


3073 






486 


836 


1128 


1566 


1953* 


2142* 


2273 


2372 


2513* 


2679* 


2828 


2966 


3078 






491 


845* 


1129* 


1588 


1968 


2145 


2288 


2377* 


2514 


2690 


2834 


2973 


3087* 


















(c) Class III (107 bursts) 














105 


493 


752 


1120 


1298 


1492 


1974 


2207 


2381 


2458 


2750 


2880 


3028 


3160 




108 


501 


753 


1125 


1306 


1634 


2035 


2230 


2382 


2460 


2760 


2917 


3068 


3166 




179 


516 


755 


1145 


1346 


1637 


2056 


2254 


2393 


2515 


2776 


2944 


3088 


3167 




228 


526 


834 


1153 


1382 


1664 


2105 


2283 


2401 


2633 


2830 


2945 


3096 






373 


555 


914 


1167 


1416 


1679 


2114 


2347 


2423 


2641 


2844 


2951 


3127 






401 


680 


942 


1190 


1435 


1693 


2129 


2349 


2424 


2677 


2848 


2980 


3139 






414 


690 


974 


1204 


1439 


1701 


2133 


2358 


2442 


2680 


2850 


2986 


3144 






465 


734t 


1114 


1221 


1443 


1747 


2152 


2368 


2453 


2719 


2861 


2990 


3146 







(d) Class IV (1 burst) 
2757* ^ 



* Placed into Class III by Gaussian model-based clustering procedure 
t Placed into Class II by Gaussian model-based clustering procedure 
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Table 5: Multivariate analysis of variance statistical tests 



Classes Wilks' Pillai's Hotelling's Wilks' A* 

A* Trace Trace F Num. dof Den. dof Prob. 

I, II, III 0.153 0.934 Im 2A5. 10 1578 <0.0001 

I, II 0.159 0.840 5.27 722. 5 683 <0.0001 

I, III 0.515 0.485 0.94 111. 5 587 <0.0001 

II, III 0.301 0.699 2.32 141. 5 304 <0.0001 



Table 6: Class properties 



Variable 


Method^ 




Class 




I 


II 


III 




(a) Means and standard 


deviations 




log T50 


NP 


1.13 ± 0.44 


-0.80 ± 0.41 


0.33 ± 0.27 


log T90 


NP 


1.55 ± 0.40 


-0.42 ± 0.44 


0.71 ± 0.32 




MB 


1.22 ± 0.39 


-0.91 ± 0.35 


0.29 ± 0.41 


log Ftot 


NP 


-5.21 ± 0.59 


-6.37 ± 0.57 


-6.11 ± 0.37 




MB 


-5.13 ± 0.58 


-6.46 ± 0.54 


-5.93 ± 0.47 


log H32I 


NP 


0.19 ± 0.27 


0.51 ± 0.27 


0.08 ± 0.40 




MB 


0.21 ± 0.26 


0.52 ± 0.28 


0.16 ± 0.40 


log Hs2 


NP 


0.43 ± 0.23 


0.70 ± 0.26 


0.35 ± 0.39 






(b) Isotropy 






<cos d> 


NP 


0.015 


-0.041 


0.010 


<sin2 6 - 1/3 > 


NP 


-0.012 


-0.025 


0.028 


Rayleigh- Watson 


NP 


0.39 


1.28 


1.80 


Bingham 


NP 


2.02 


7.32 


1.95 






(c) Summary 






Number 


NP 


486 


203 


107 




MB 


426 


170 


201 


Duration 




long 


short 


intermediate 


Fluence 




bright 


faint 


intermediate 


Spectrum 




intermediate 


hard 


soft 



NP = nonparametric clustering analysis in 5 dimensions (§4) 
MB = model-based clustering analysis in 3 dimensions (§5) 
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Table 7: Correlation coefficients within classes 



log T50 log T90 log Ftot log F321 log H32 



(a) Class I - 


Nonparametric clustering (N=486) 


|r| > 0. 


15 corresponds to P < 


0.001 significance level 


log T50 


1.00 






log Tgo 


0.88 


1.00 




log Ftot 


0.10 


0.22 


1.00 


log i?321 


-0.11 


-0.08 


0.39 1.00 


log -ff32 


-0.11 


-0.08 


0.38 0.97 1.00 


(b) 


Class I - 


- Model-based 


clustering (N=426) 


\r\ >0. 


16 corresponds to P < 


0.001 significance level 


log Tso 


1.00 






log Tgo 


N/A 


N/A 




log Ftot 


N/A 


0.01 


1.00 


log if321 


N/A 


-0.01 


0.06 1.00 


log -^32 


N/A 


N/A 


N/A N/A N/A 


(d) Class II 


- Model-based clustering (N=170) 


r > 0.25 corresponds to P < 0.001 significance level 


log T50 


1.00 






log Tgo 


N/A 


N/A 




log Ftot 


N/A 


-0.03 


1.00 


log H321 


N/A 


-0.05 


0.02 1.00 


log Hs2 


N/A 


N/A 


N/A N/A N/A 


(e) Class III - 


- Nonparametric clustering (N=107) 


|r| > 0.32 corrc; 


sponds to P < 0.001 significance level 


log Tso 


1.00 






log Tgo 


0.86 


1.00 




log Ftot 


0.02 


0.06 


1.00 


log H32I 


-0.24 


-0.34 


-0.16 1.00 


log H32 


-0.22 


-0.32 


-0.22 0.95 1.00 


(f) Class III 


- Model-based clustering (N=201) 


r > 0.23 corresponds to P < 0.001 significance level 


log Tso 


1.00 






log T90 


N/A 


N/A 




log Ftot 


N/A 


0.03 


1.00 


log ^321 


N/A 


-0.01 


0.07 1.00 


log H32 


N/A 


N/A 


N/A N/A N/A 



